FREC 5984 Project: Fit to Historical Data
VT NEE Team: Josh Rady, Marcus Taylor, and Tyler Weiglein
01 Apr 2021
Read in, format, and visualize data.
# Read in target data
target_data <- read.csv("terrestrial_30min-targets.csv")
# Read in gap-filled temperature data
temp_data <- read.csv("temp_gapfill.csv")
# Convert time column from character to POSIXct
# Time zone in target data is UTC according to https://docs.google.com/document/d/1l7sxBk-z-GHTlk50rdxP0lPTwJzFJ2gykclINkMsWcc/edit
target_data$time <- ymd_hms(target_data$time, tz = "UTC")
temp_data$startDateTime <- ymd_hms(temp_data$startDateTime, tz = "UTC") # Ask Marcus to verify UTC
# Convert siteID column from character to factor
target_data$siteID <- as.factor(target_data$siteID)
temp_data$siteID <- as.factor(temp_data$siteID)
# Plot all data
ggplot(data = target_data) +
geom_point(aes(x = time, y = nee, color = siteID)) +
facet_grid(siteID ~ .) +
labs(title = "Full Time Period: 01 Feb 2017 - 31 Jan 2021",
x = "Time (UTC)",
y = expression(paste("NEE (",mu, "mol" ~ CO[2] ~ m^-2 ~ s^-1,")"))) +
theme_bw() +
theme(legend.position = "none")
# Plot data for example week (2020-06-07 through 2020-06-13)
ggplot(data = filter(target_data,
time >= "2020-06-07",
time < "2020-06-14")) +
geom_point(aes(x = time, y = nee, color = siteID)) +
facet_grid(siteID ~ .) +
labs(title = "Example Time Period: 07 Jun 2020 - 13 Jun 2020",
x = "Time (UTC)",
y = expression(paste("NEE (",mu, "mol" ~ CO[2] ~ m^-2 ~ s^-1,")"))) +
theme_bw() +
theme(legend.position = "none")
Construct null model (random walk) for NIMBLE.
RandomWalk <- nimbleCode({
# Note:
# x = real NEE (state variable)
# y = observed NEE
# Priors
x[1] ~ dnorm(x_ic, sd = sd_ic)
sd_add ~ dunif(0, 100)
# Process model
for(t in 2:n){
pred[t] <- x[t-1]
x[t] ~ dnorm(pred[t], sd = sd_add)
}
# Data model
for(t in 1:n){
y[t] ~ dnorm(x[t], sd = sd_obs[t])
}
})Construct temperature dynamic linear model (DLM) for NIMBLE.
TempDLM <- nimbleCode({
# Note:
# x = real NEE (state variable)
# y = observed NEE
# Priors
x[1] ~ dnorm(x_ic, sd = sd_ic)
sd_add ~ dunif(0, 100)
b1 ~ dnorm(0, 100) # Temperature effect coefficient
# Process model
for(t in 2:n){
pred[t] <- x[t-1] + b1 * temp[t]
x[t] ~ dnorm(pred[t], sd = sd_add)
}
# Data model
for(t in 1:n){
y[t] ~ dnorm(x[t], sd = sd_obs[t])
}
})Because of the size of the data set, a subset of the data had to be used for model fitting. The variables controlling the size of the data subset are specified below.
# Note: B/c of size of data set (and potentially number of NAs), need to use subset of data
subset_length_days <- 7 # Length of subset [d]
subset_length_timesteps <- subset_length_days * 24 * 2 # Number of time steps in subset (days x hours/day x half hours/hour)Specify data from BART to be used for model fitting. (Because of the large difference between NEE values in the winter versus other times of year, the additional restriction of limiting the data subset to April-November was used.)
# Separate data by site
target_data_BART <- filter(target_data, siteID == "BART")
temp_data_BART <- filter(temp_data, siteID == "BART")
# Find data subset w/ fewest number of NAs
min_n_NAs <- subset_length_timesteps # Initialize scalar for fewest NAs
BART_start_i <- NA # Initialize BART starting index
# Loop through data to find data subset (not during winter) with minimum number of NAs for specified subset length
for(i in 1:(length(target_data_BART$nee) - subset_length_timesteps + 1)){
num_NAs <- sum(is.na(target_data_BART$nee[i:(i + subset_length_timesteps - 1)])) # Count number of NEE NAs in data subset
if(num_NAs <= min_n_NAs & sum(substr(target_data_BART$time[i], 6, 7) != c("12", "01", "02", "03")) == 4){ # Note: Using <= so that most recent data is used; also, excluding data from December through March
BART_start_i <- i # Save new starting index
min_n_NAs <- num_NAs # Save new lowest number of NAs
}
}
# Calculate data subset end index
BART_end_i <- BART_start_i + subset_length_timesteps - 1
# Subset data
target_data_BART_sub <- target_data_BART[BART_start_i:BART_end_i,]
# Interpolate NEE values to estimate sd_obs
nee_interp_BART_sub <- approx(x = target_data_BART_sub$time[!is.na(target_data_BART_sub$nee)],
y = target_data_BART_sub$nee[!is.na(target_data_BART_sub$nee)],
xout = target_data_BART_sub$time,
method = "linear",
yleft = mean(target_data_BART_sub$nee, na.rm = TRUE),
yright = mean(target_data_BART_sub$nee, na.rm = TRUE),
na.rm = TRUE)
sd_obs_BART_sub <- rep(NA, length(target_data_BART_sub$nee)) # Initialize sd_obs vector
for(i in 1:length(sd_obs_BART_sub)){
if(nee_interp_BART_sub$y[i] >= 0){ # Calculate sd_obs if NEE is positive
sd_obs_BART_sub[i] <- target_data_BART_sub$nee_sd_intercept[1] + target_data_BART_sub$ nee_sd_slopeP[1] * nee_interp_BART_sub$y[i] # Note: Since intercept and slope do not change, can just use first value of intercept and slope values
}else{ # Calculate sd_obs if NEE is negative
sd_obs_BART_sub[i] <- target_data_BART_sub$nee_sd_intercept[1] + target_data_BART_sub$ nee_sd_slopeN[1] * nee_interp_BART_sub$y[i]
}
}
# Match temperature data to subset NEE data
BART_temp_i <- match(target_data_BART_sub$time, temp_data_BART$startDateTime)
temp_data_BART_sub <- temp_data_BART[BART_temp_i,]
# Specify data (Note: RW = random walk)
data_RW_BART <- list(y = target_data_BART_sub$nee,
sd_obs = sd_obs_BART_sub)
data_TempDLM_BART <- list(y = target_data_BART_sub$nee,
temp = temp_data_BART_sub$temp,
sd_obs = sd_obs_BART_sub)Specify data from KONZ to be used for model fitting.
# Separate data by site
target_data_KONZ <- filter(target_data, siteID == "KONZ")
temp_data_KONZ <- filter(temp_data, siteID == "KONZ")
# Find data subset w/ fewest number of NAs
min_n_NAs <- subset_length_timesteps # Initialize scalar for fewest NAs
KONZ_start_i <- NA # Initialize KONZ starting index
# Loop through data to find data subset with minimum number of NAs for specified subset length
for(i in 1:(length(target_data_KONZ$nee) - subset_length_timesteps + 1)){
num_NAs <- sum(is.na(target_data_KONZ$nee[i:(i + subset_length_timesteps - 1)])) # Count number of NEE NAs in data subset
if(num_NAs <= min_n_NAs){ # Note: Using <= so that most recent data is used
KONZ_start_i <- i # Save new starting index
min_n_NAs <- num_NAs # Save new lowest number of NAs
}
}
# Calculate data subset end index
KONZ_end_i <- KONZ_start_i + subset_length_timesteps - 1
# Subset data
target_data_KONZ_sub <- target_data_KONZ[KONZ_start_i:KONZ_end_i,]
# Interpolate NEE values to estimate sd_obs
nee_interp_KONZ_sub <- approx(x = target_data_KONZ_sub$time[!is.na(target_data_KONZ_sub$nee)],
y = target_data_KONZ_sub$nee[!is.na(target_data_KONZ_sub$nee)],
xout = target_data_KONZ_sub$time,
method = "linear",
yleft = mean(target_data_KONZ_sub$nee, na.rm = TRUE),
yright = mean(target_data_KONZ_sub$nee, na.rm = TRUE),
na.rm = TRUE)
sd_obs_KONZ_sub <- rep(NA, length(target_data_KONZ_sub$nee)) # Initialize sd_obs vector
for(i in 1:length(sd_obs_KONZ_sub)){
if(nee_interp_KONZ_sub$y[i] >= 0){ # Calculate sd_obs if NEE is positive
sd_obs_KONZ_sub[i] <- target_data_KONZ_sub$nee_sd_intercept[1] + target_data_KONZ_sub$ nee_sd_slopeP[1] * nee_interp_KONZ_sub$y[i] # Note: Since intercept and slope do not change, can just use first value of intercept and slope values
}else{ # Calculate sd_obs if NEE is negative
sd_obs_KONZ_sub[i] <- target_data_KONZ_sub$nee_sd_intercept[1] + target_data_KONZ_sub$ nee_sd_slopeN[1] * nee_interp_KONZ_sub$y[i]
}
}
# Match temperature data to subset NEE data
KONZ_temp_i <- match(target_data_KONZ_sub$time, temp_data_KONZ$startDateTime)
temp_data_KONZ_sub <- temp_data_KONZ[KONZ_temp_i,]
# Specify data (Note: RW = random walk)
data_RW_KONZ <- list(y = target_data_KONZ_sub$nee,
sd_obs = sd_obs_KONZ_sub)
data_TempDLM_KONZ <- list(y = target_data_KONZ_sub$nee,
temp = temp_data_KONZ_sub$temp,
sd_obs = sd_obs_KONZ_sub)Specify data from OSBS to be used for model fitting.
# Separate data by site
target_data_OSBS <- filter(target_data, siteID == "OSBS")
temp_data_OSBS <- filter(temp_data, siteID == "OSBS")
# Remove duplicate rows
target_data_OSBS <- unique(target_data_OSBS)
# Find data subset w/ fewest number of NAs
min_n_NAs <- subset_length_timesteps # Initialize scalar for fewest NAs
OSBS_start_i <- NA # Initialize OSBS starting index
# Loop through data to find data subset with minimum number of NAs for specified subset length
for(i in 1:(length(target_data_OSBS$nee) - subset_length_timesteps + 1)){
num_NAs <- sum(is.na(target_data_OSBS$nee[i:(i + subset_length_timesteps - 1)])) # Count number of NEE NAs in data subset
if(num_NAs <= min_n_NAs){ # Note: Using <= so that most recent data is used
OSBS_start_i <- i # Save new starting index
min_n_NAs <- num_NAs # Save new lowest number of NAs
}
}
# Calculate data subset end index
OSBS_end_i <- OSBS_start_i + subset_length_timesteps - 1
# Subset data
target_data_OSBS_sub <- target_data_OSBS[OSBS_start_i:OSBS_end_i,]
# Interpolate NEE values to estimate sd_obs
nee_interp_OSBS_sub <- approx(x = target_data_OSBS_sub$time[!is.na(target_data_OSBS_sub$nee)],
y = target_data_OSBS_sub$nee[!is.na(target_data_OSBS_sub$nee)],
xout = target_data_OSBS_sub$time,
method = "linear",
yleft = mean(target_data_OSBS_sub$nee, na.rm = TRUE),
yright = mean(target_data_OSBS_sub$nee, na.rm = TRUE),
na.rm = TRUE)
sd_obs_OSBS_sub <- rep(NA, length(target_data_OSBS_sub$nee)) # Initialize sd_obs vector
for(i in 1:length(sd_obs_OSBS_sub)){
if(nee_interp_OSBS_sub$y[i] >= 0){ # Calculate sd_obs if NEE is positive
sd_obs_OSBS_sub[i] <- target_data_OSBS_sub$nee_sd_intercept[1] + target_data_OSBS_sub$ nee_sd_slopeP[1] * nee_interp_OSBS_sub$y[i] # Note: Since intercept and slope do not change, can just use first value of intercept and slope values
}else{ # Calculate sd_obs if NEE is negative
sd_obs_OSBS_sub[i] <- target_data_OSBS_sub$nee_sd_intercept[1] + target_data_OSBS_sub$ nee_sd_slopeN[1] * nee_interp_OSBS_sub$y[i]
}
}
# Match temperature data to subset NEE data
OSBS_temp_i <- match(target_data_OSBS_sub$time, temp_data_OSBS$startDateTime)
temp_data_OSBS_sub <- temp_data_OSBS[OSBS_temp_i,]
# Specify data (Note: RW = random walk)
data_RW_OSBS <- list(y = target_data_OSBS_sub$nee,
sd_obs = sd_obs_OSBS_sub)
data_TempDLM_OSBS <- list(y = target_data_OSBS_sub$nee,
temp = temp_data_OSBS_sub$temp,
sd_obs = sd_obs_OSBS_sub)Specify data from SRER to be used for model fitting.
# Separate data by site
target_data_SRER <- filter(target_data, siteID == "SRER")
temp_data_SRER <- filter(temp_data, siteID == "SRER")
# Remove duplicate rows
target_data_SRER <- unique(target_data_SRER)
# Find data subset w/ fewest number of NAs
min_n_NAs <- subset_length_timesteps # Initialize scalar for fewest NAs
SRER_start_i <- NA # Initialize SRER starting index
# Loop through data to find data subset with minimum number of NAs for specified subset length
for(i in 1:(length(target_data_SRER$nee) - subset_length_timesteps + 1)){
num_NAs <- sum(is.na(target_data_SRER$nee[i:(i + subset_length_timesteps - 1)])) # Count number of NEE NAs in data subset
if(num_NAs <= min_n_NAs){ # Note: Using <= so that most recent data is used
SRER_start_i <- i # Save new starting index
min_n_NAs <- num_NAs # Save new lowest number of NAs
}
}
# Calculate data subset end index
SRER_end_i <- SRER_start_i + subset_length_timesteps - 1
# Subset data
target_data_SRER_sub <- target_data_SRER[SRER_start_i:SRER_end_i,]
# Interpolate NEE values to estimate sd_obs
nee_interp_SRER_sub <- approx(x = target_data_SRER_sub$time[!is.na(target_data_SRER_sub$nee)],
y = target_data_SRER_sub$nee[!is.na(target_data_SRER_sub$nee)],
xout = target_data_SRER_sub$time,
method = "linear",
yleft = mean(target_data_SRER_sub$nee, na.rm = TRUE),
yright = mean(target_data_SRER_sub$nee, na.rm = TRUE),
na.rm = TRUE)
sd_obs_SRER_sub <- rep(NA, length(target_data_SRER_sub$nee)) # Initialize sd_obs vector
for(i in 1:length(sd_obs_SRER_sub)){
if(nee_interp_SRER_sub$y[i] >= 0){ # Calculate sd_obs if NEE is positive
sd_obs_SRER_sub[i] <- target_data_SRER_sub$nee_sd_intercept[1] + target_data_SRER_sub$ nee_sd_slopeP[1] * nee_interp_SRER_sub$y[i] # Note: Since intercept and slope do not change, can just use first value of intercept and slope values
}else{ # Calculate sd_obs if NEE is negative
sd_obs_SRER_sub[i] <- target_data_SRER_sub$nee_sd_intercept[1] + target_data_SRER_sub$ nee_sd_slopeN[1] * nee_interp_SRER_sub$y[i]
}
}
# Match temperature data to subset NEE data
SRER_temp_i <- match(target_data_SRER_sub$time, temp_data_SRER$startDateTime)
temp_data_SRER_sub <- temp_data_SRER[SRER_temp_i,]
# Specify data (Note: RW = random walk)
data_RW_SRER <- list(y = target_data_SRER_sub$nee,
sd_obs = sd_obs_SRER_sub)
data_TempDLM_SRER <- list(y = target_data_SRER_sub$nee,
temp = temp_data_SRER_sub$temp,
sd_obs = sd_obs_SRER_sub)Specify constants for NIMBLE.
# Specify constants
constants_BART <- list(n = length(target_data_BART_sub$nee),
x_ic = 0,
sd_ic = 1)
constants_KONZ <- list(n = length(target_data_KONZ_sub$nee),
x_ic = 0,
sd_ic = 1)
constants_OSBS <- list(n = length(target_data_OSBS_sub$nee),
x_ic = 0,
sd_ic = 1)
constants_SRER <- list(n = length(target_data_SRER_sub$nee),
x_ic = 0,
sd_ic = 1)Specify initial conditions for null (random walk) model.
# Specify number of chains
nchains = 3
# Initialize initial condition lists (Note: RW = random walk)
inits_RW_BART <- list()
inits_RW_KONZ <- list()
inits_RW_OSBS <- list()
inits_RW_SRER <- list()
# Generate initial conditions
for(i in 1:nchains){
# BART
y.samp <- sample(nee_interp_BART_sub$y, length(nee_interp_BART_sub$y), replace = TRUE)
inits_RW_BART[[i]] <- list(sd_add = sd(diff(y.samp), na.rm = TRUE),
x = nee_interp_BART_sub$y)
# KONZ
y.samp <- sample(nee_interp_KONZ_sub$y, length(nee_interp_KONZ_sub$y), replace = TRUE)
inits_RW_KONZ[[i]] <- list(sd_add = sd(diff(y.samp), na.rm = TRUE),
x = nee_interp_KONZ_sub$y)
# OSBS
y.samp <- sample(nee_interp_OSBS_sub$y, length(nee_interp_OSBS_sub$y), replace = TRUE)
inits_RW_OSBS[[i]] <- list(sd_add = sd(diff(y.samp), na.rm = TRUE),
x = nee_interp_OSBS_sub$y)
# SRER
y.samp <- sample(nee_interp_SRER_sub$y, length(nee_interp_SRER_sub$y), replace = TRUE)
inits_RW_SRER[[i]] <- list(sd_add = sd(diff(y.samp), na.rm = TRUE),
x = nee_interp_SRER_sub$y)
}Specify initial conditions for temperature DLM.
# Initialize initial condition lists
inits_TempDLM_BART <- list()
inits_TempDLM_KONZ <- list()
inits_TempDLM_OSBS <- list()
inits_TempDLM_SRER <- list()
# Generate initial conditions
for(i in 1:nchains){
# BART
y.samp <- sample(nee_interp_BART_sub$y, length(nee_interp_BART_sub$y), replace = TRUE)
inits_TempDLM_BART[[i]] <- list(sd_add = sd(diff(y.samp), na.rm = TRUE),
x = nee_interp_BART_sub$y,
b1 = rnorm(1, mean = 0, sd = 1))
# KONZ
y.samp <- sample(nee_interp_KONZ_sub$y, length(nee_interp_KONZ_sub$y), replace = TRUE)
inits_TempDLM_KONZ[[i]] <- list(sd_add = sd(diff(y.samp), na.rm = TRUE),
x = nee_interp_KONZ_sub$y,
b1 = rnorm(1, mean = 0, sd = 1))
# OSBS
y.samp <- sample(nee_interp_OSBS_sub$y, length(nee_interp_OSBS_sub$y), replace = TRUE)
inits_TempDLM_OSBS[[i]] <- list(sd_add = sd(diff(y.samp), na.rm = TRUE),
x = nee_interp_OSBS_sub$y,
b1 = rnorm(1, mean = 0, sd = 1))
# SRER
y.samp <- sample(nee_interp_SRER_sub$y, length(nee_interp_SRER_sub$y), replace = TRUE)
inits_TempDLM_SRER[[i]] <- list(sd_add = sd(diff(y.samp), na.rm = TRUE),
x = nee_interp_SRER_sub$y,
b1 = rnorm(1, mean = 0, sd = 1))
}Run NIMBLE for null (random walk) model.
# BART
nimble_out_RW_BART <- nimbleMCMC(code = RandomWalk,
data = data_RW_BART,
inits = inits_RW_BART,
constants = constants_BART,
monitors = c("sd_add",
"x",
"y"),
niter = 11000,
nchains = nchains,
samplesAsCodaMCMC = TRUE)## NAs were detected in model variables: y, logProb_y.
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|# KONZ
nimble_out_RW_KONZ <- nimbleMCMC(code = RandomWalk,
data = data_RW_KONZ,
inits = inits_RW_KONZ,
constants = constants_KONZ,
monitors = c("sd_add",
"x",
"y"),
niter = 11000,
nchains = nchains,
samplesAsCodaMCMC = TRUE)## NAs were detected in model variables: y, logProb_y.
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|# OSBS
nimble_out_RW_OSBS <- nimbleMCMC(code = RandomWalk,
data = data_RW_OSBS,
inits = inits_RW_OSBS,
constants = constants_OSBS,
monitors = c("sd_add",
"x",
"y"),
niter = 11000,
nchains = nchains,
samplesAsCodaMCMC = TRUE)## NAs were detected in model variables: y, logProb_y.
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|# SRER
nimble_out_RW_SRER <- nimbleMCMC(code = RandomWalk,
data = data_RW_SRER,
inits = inits_RW_SRER,
constants = constants_SRER,
monitors = c("sd_add",
"x",
"y"),
niter = 11000,
nchains = nchains,
samplesAsCodaMCMC = TRUE)## NAs were detected in model variables: y, logProb_y.
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|Run NIMBLE for temperature DLM.
# BART
nimble_out_TempDLM_BART <- nimbleMCMC(code = TempDLM,
data = data_TempDLM_BART,
inits = inits_TempDLM_BART,
constants = constants_BART,
monitors = c("b1",
"sd_add",
"x",
"y"),
niter = 11000,
nchains = nchains,
samplesAsCodaMCMC = TRUE)## NAs were detected in model variables: y, logProb_y.
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|# KONZ
nimble_out_TempDLM_KONZ <- nimbleMCMC(code = TempDLM,
data = data_TempDLM_KONZ,
inits = inits_TempDLM_KONZ,
constants = constants_KONZ,
monitors = c("b1",
"sd_add",
"x",
"y"),
niter = 11000,
nchains = nchains,
samplesAsCodaMCMC = TRUE)## NAs were detected in model variables: y, logProb_y.
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|# OSBS
nimble_out_TempDLM_OSBS <- nimbleMCMC(code = TempDLM,
data = data_TempDLM_OSBS,
inits = inits_TempDLM_OSBS,
constants = constants_OSBS,
monitors = c("b1",
"sd_add",
"x",
"y"),
niter = 11000,
nchains = nchains,
samplesAsCodaMCMC = TRUE)## NAs were detected in model variables: y, logProb_y.
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|# SRER
nimble_out_TempDLM_SRER <- nimbleMCMC(code = TempDLM,
data = data_TempDLM_SRER,
inits = inits_TempDLM_SRER,
constants = constants_SRER,
monitors = c("b1",
"sd_add",
"x",
"y"),
niter = 11000,
nchains = nchains,
samplesAsCodaMCMC = TRUE)## NAs were detected in model variables: y, logProb_y.
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|Visualize chains and remove burn-in for null (random walk) model.
# BART
plot(nimble_out_RW_BART[, "sd_add"], main = "BART: sd_add") # Visualize all chains
burnin_RW_BART <- 2000 # Burn-in
nimble_burn_RW_BART <- window(nimble_out_RW_BART, start = burnin_RW_BART) # Exclude burn-in
plot(nimble_burn_RW_BART[, "sd_add"], main = "BART: sd_add (burn-in removed)") # Visualize chains w/o burn-in
# KONZ
plot(nimble_out_RW_KONZ[, "sd_add"], main = "KONZ: sd_add") # Visualize all chains
burnin_RW_KONZ <- 2000 # Burn-in
nimble_burn_RW_KONZ <- window(nimble_out_RW_KONZ, start = burnin_RW_KONZ) # Exclude burn-in
plot(nimble_burn_RW_KONZ[, "sd_add"], main = "KONZ: sd_add (burn-in removed)") # Visualize chains w/o burn-in
# OSBS
plot(nimble_out_RW_OSBS[, "sd_add"], main = "OSBS: sd_add") # Visualize all chains
burnin_RW_OSBS <- 2000 # Burn-in
nimble_burn_RW_OSBS <- window(nimble_out_RW_OSBS, start = burnin_RW_OSBS) # Exclude burn-in
plot(nimble_burn_RW_OSBS[, "sd_add"], main = "OSBS: sd_add (burn-in removed)") # Visualize chains w/o burn-in
# SRER
plot(nimble_out_RW_SRER[, "sd_add"], main = "SRER: sd_add") # Visualize all chains
burnin_RW_SRER <- 2000 # Burn-in
nimble_burn_RW_SRER <- window(nimble_out_RW_SRER, start = burnin_RW_SRER) # Exclude burn-in
plot(nimble_burn_RW_SRER[, "sd_add"], main = "SRER: sd_add (burn-in removed)") # Visualize chains w/o burn-in
Visualize chains and remove burn-in for temperature DLM.
# BART
plot(nimble_out_TempDLM_BART[, c("b1")], main = "BART: b1") # Visualize all chains
plot(nimble_out_TempDLM_BART[, c("sd_add")], main = "BART: sd_add") # Visualize all chains
burnin_TempDLM_BART <- 2000 # Burn-in
nimble_burn_TempDLM_BART <- window(nimble_out_TempDLM_BART, start = burnin_TempDLM_BART) # Exclude burn-in
plot(nimble_burn_TempDLM_BART[, c("b1")], main = "BART: b1 (burn-in removed)") # Visualize chains w/o burn-in
plot(nimble_burn_TempDLM_BART[, c("sd_add")], main = "BART: sd_add (burn-in removed)") # Visualize chains w/o burn-in
# KONZ
plot(nimble_out_TempDLM_KONZ[, c("b1")], main = "KONZ: b1") # Visualize all chains
plot(nimble_out_TempDLM_KONZ[, c("sd_add")], main = "KONZ: sd_add") # Visualize all chains
burnin_TempDLM_KONZ <- 2000 # Burn-in
nimble_burn_TempDLM_KONZ <- window(nimble_out_TempDLM_KONZ, start = burnin_TempDLM_KONZ) # Exclude burn-in
plot(nimble_burn_TempDLM_KONZ[, c("b1")], main = "KONZ: b1 (burn-in removed)") # Visualize chains w/o burn-in
plot(nimble_burn_TempDLM_KONZ[, c("sd_add")], main = "KONZ: sd_add (burn-in removed)") # Visualize chains w/o burn-in
# OSBS
plot(nimble_out_TempDLM_OSBS[, c("b1")], main = "OSBS: b1") # Visualize all chains
plot(nimble_out_TempDLM_OSBS[, c("sd_add")], main = "OSBS: sd_add") # Visualize all chains
burnin_TempDLM_OSBS <- 2000 # Burn-in
nimble_burn_TempDLM_OSBS <- window(nimble_out_TempDLM_OSBS, start = burnin_TempDLM_OSBS) # Exclude burn-in
plot(nimble_burn_TempDLM_OSBS[, c("b1")], main = "OSBS: b1 (burn-in removed)") # Visualize chains w/o burn-in
plot(nimble_burn_TempDLM_OSBS[, c("sd_add")], main = "OSBS: sd_add (burn-in removed)") # Visualize chains w/o burn-in
# SRER
plot(nimble_out_TempDLM_SRER[, c("b1")], main = "SRER: b1") # Visualize all chains
plot(nimble_out_TempDLM_SRER[, c("sd_add")], main = "SRER: sd_add") # Visualize all chains
burnin_TempDLM_SRER <- 2000 # Burn-in
nimble_burn_TempDLM_SRER <- window(nimble_out_TempDLM_SRER, start = burnin_TempDLM_SRER) # Exclude burn-in
plot(nimble_burn_TempDLM_SRER[, c("b1")], main = "SRER: b1 (burn-in removed)") # Visualize chains w/o burn-in
plot(nimble_burn_TempDLM_SRER[, c("sd_add")], main = "SRER: sd_add (burn-in removed)") # Visualize chains w/o burn-in
Examine diagnostics for null (random walk) model. Since the upper confidence interval for the Brooks-Gelman-Rubin convergence diagnostic is less than 1.1, we can conclude the chains have converged.
# BART
gelman.diag(nimble_burn_RW_BART[, "sd_add"])## Potential scale reduction factors:
##
## Point est. Upper C.I.
## [1,] 1 1.01# KONZ
gelman.diag(nimble_burn_RW_KONZ[, "sd_add"])## Potential scale reduction factors:
##
## Point est. Upper C.I.
## [1,] 1 1.01# OSBS
gelman.diag(nimble_burn_RW_OSBS[, "sd_add"])## Potential scale reduction factors:
##
## Point est. Upper C.I.
## [1,] 1.01 1.03# SRER
gelman.diag(nimble_burn_RW_SRER[, "sd_add"])## Potential scale reduction factors:
##
## Point est. Upper C.I.
## [1,] 1.01 1.03Examine diagnostics for temperature DLM. Since the upper confidence interval for the Brooks-Gelman-Rubin convergence diagnostic is less than 1.1, we can conclude the chains have converged.
# BART
gelman.diag(nimble_burn_TempDLM_BART[, "b1"])## Potential scale reduction factors:
##
## Point est. Upper C.I.
## [1,] 1.01 1.02gelman.diag(nimble_burn_TempDLM_BART[, "sd_add"])## Potential scale reduction factors:
##
## Point est. Upper C.I.
## [1,] 1.01 1.03# KONZ
gelman.diag(nimble_burn_TempDLM_KONZ[, "b1"])## Potential scale reduction factors:
##
## Point est. Upper C.I.
## [1,] 1 1gelman.diag(nimble_burn_TempDLM_KONZ[, "sd_add"])## Potential scale reduction factors:
##
## Point est. Upper C.I.
## [1,] 1 1.01# OSBS
gelman.diag(nimble_burn_TempDLM_OSBS[, "b1"])## Potential scale reduction factors:
##
## Point est. Upper C.I.
## [1,] 1 1gelman.diag(nimble_burn_TempDLM_OSBS[, "sd_add"])## Potential scale reduction factors:
##
## Point est. Upper C.I.
## [1,] 1.01 1.04# SRER
gelman.diag(nimble_burn_TempDLM_SRER[, "b1"])## Potential scale reduction factors:
##
## Point est. Upper C.I.
## [1,] 1 1gelman.diag(nimble_burn_TempDLM_SRER[, "sd_add"])## Potential scale reduction factors:
##
## Point est. Upper C.I.
## [1,] 1 1.01Convert null (random walk) model output using tidybayes.
# BART
chain_RW_BART <- nimble_burn_RW_BART %>%
spread_draws(y[day], x[day], sd_add)
# KONZ
chain_RW_KONZ <- nimble_burn_RW_KONZ %>%
spread_draws(y[day], x[day], sd_add)
# OSBS
chain_RW_OSBS <- nimble_burn_RW_OSBS %>%
spread_draws(y[day], x[day], sd_add)
# SRER
chain_RW_SRER <- nimble_burn_RW_SRER %>%
spread_draws(y[day], x[day], sd_add)Convert temperature DLM output using tidybayes.
# BART
chain_TempDLM_BART <- nimble_burn_TempDLM_BART %>%
spread_draws(y[day], x[day], sd_add)
# KONZ
chain_TempDLM_KONZ <- nimble_burn_TempDLM_KONZ %>%
spread_draws(y[day], x[day], sd_add)
# OSBS
chain_TempDLM_OSBS <- nimble_burn_TempDLM_OSBS %>%
spread_draws(y[day], x[day], sd_add)
# SRER
chain_TempDLM_SRER <- nimble_burn_TempDLM_SRER %>%
spread_draws(y[day], x[day], sd_add)Plot null (random walk) model versus observations.
plot_colors <- c("#F8766D", "#7CAE00", "#00BFC4", "#C77CFF")
# BART
plot_chain_RW_BART <- chain_RW_BART %>%
group_by(day) %>%
summarise(mean = mean(x),
lower = quantile(x, 0.025),
upper = quantile(x, 0.975),
.groups = "drop") %>%
mutate(time = target_data_BART_sub$time)
plot_RW_BART <- ggplot(data = plot_chain_RW_BART, aes(x = time, y = mean)) +
geom_line() +
geom_ribbon(aes(ymin = lower, ymax = upper),
alpha = 0.2,
color = plot_colors[1],
fill = plot_colors[1]) +
geom_point(data = target_data_BART_sub, aes(x = time, y = nee), color = "black") +
labs(title = paste("Site: BART | Time Period: ", substr(target_data_BART_sub$time[1], 1, 10), "to", substr(target_data_BART_sub$time[subset_length_timesteps], 1, 10)),
x = "Time (UTC)",
y = expression(paste("NEE (",mu, "mol" ~ CO[2] ~ m^-2 ~ s^-1,")"))) +
theme_bw()
# KONZ
plot_chain_RW_KONZ <- chain_RW_KONZ %>%
group_by(day) %>%
summarise(mean = mean(x),
lower = quantile(x, 0.025),
upper = quantile(x, 0.975),
.groups = "drop") %>%
mutate(time = target_data_KONZ_sub$time)
plot_RW_KONZ <- ggplot(data = plot_chain_RW_KONZ, aes(x = time, y = mean)) +
geom_line() +
geom_ribbon(aes(ymin = lower, ymax = upper),
alpha = 0.2,
color = plot_colors[2],
fill = plot_colors[2]) +
geom_point(data = target_data_KONZ_sub, aes(x = time, y = nee), color = "black") +
labs(title = paste("Site: KONZ | Time Period: ", substr(target_data_KONZ_sub$time[1], 1, 10), "to", substr(target_data_KONZ_sub$time[subset_length_timesteps], 1, 10)),
x = "Time (UTC)",
y = expression(paste("NEE (",mu, "mol" ~ CO[2] ~ m^-2 ~ s^-1,")"))) +
theme_bw()
# OSBS
plot_chain_RW_OSBS <- chain_RW_OSBS %>%
group_by(day) %>%
summarise(mean = mean(x),
lower = quantile(x, 0.025),
upper = quantile(x, 0.975),
.groups = "drop") %>%
mutate(time = target_data_OSBS_sub$time)
plot_RW_OSBS <- ggplot(data = plot_chain_RW_OSBS, aes(x = time, y = mean)) +
geom_line() +
geom_ribbon(aes(ymin = lower, ymax = upper),
alpha = 0.2,
color = plot_colors[3],
fill = plot_colors[3]) +
geom_point(data = target_data_OSBS_sub, aes(x = time, y = nee), color = "black") +
labs(title = paste("Site: OSBS | Time Period: ", substr(target_data_OSBS_sub$time[1], 1, 10), "to", substr(target_data_OSBS_sub$time[subset_length_timesteps], 1, 10)),
x = "Time (UTC)",
y = expression(paste("NEE (",mu, "mol" ~ CO[2] ~ m^-2 ~ s^-1,")"))) +
theme_bw()
# SRER
plot_chain_RW_SRER <- chain_RW_SRER %>%
group_by(day) %>%
summarise(mean = mean(x),
lower = quantile(x, 0.025),
upper = quantile(x, 0.975),
.groups = "drop") %>%
mutate(time = target_data_SRER_sub$time)
plot_RW_SRER <- ggplot(data = plot_chain_RW_SRER, aes(x = time, y = mean)) +
geom_line() +
geom_ribbon(aes(ymin = lower, ymax = upper),
alpha = 0.2,
color = plot_colors[4],
fill = plot_colors[4]) +
geom_point(data = target_data_SRER_sub, aes(x = time, y = nee), color = "black") +
labs(title = paste("Site: SRER | Time Period: ", substr(target_data_SRER_sub$time[1], 1, 10), "to", substr(target_data_SRER_sub$time[subset_length_timesteps], 1, 10)),
x = "Time (UTC)",
y = expression(paste("NEE (",mu, "mol" ~ CO[2] ~ m^-2 ~ s^-1,")"))) +
theme_bw()
plot_RW_BART
plot_RW_KONZ
plot_RW_OSBS
plot_RW_SRER
Plot temperature DLM versus observations.
# BART
plot_chain_TempDLM_BART <- chain_TempDLM_BART %>%
group_by(day) %>%
summarise(mean = mean(x),
lower = quantile(x, 0.025),
upper = quantile(x, 0.975),
.groups = "drop") %>%
mutate(time = target_data_BART_sub$time)
plot_TempDLM_BART <- ggplot(data = plot_chain_TempDLM_BART, aes(x = time, y = mean)) +
geom_line() +
geom_ribbon(aes(ymin = lower, ymax = upper),
alpha = 0.2,
color = plot_colors[1],
fill = plot_colors[1]) +
geom_point(data = target_data_BART_sub, aes(x = time, y = nee), color = "black") +
labs(title = paste("Site: BART | Time Period: ", substr(target_data_BART_sub$time[1], 1, 10), "to", substr(target_data_BART_sub$time[subset_length_timesteps], 1, 10)),
x = "Time (UTC)",
y = expression(paste("NEE (",mu, "mol" ~ CO[2] ~ m^-2 ~ s^-1,")"))) +
theme_bw()
# KONZ
plot_chain_TempDLM_KONZ <- chain_TempDLM_KONZ %>%
group_by(day) %>%
summarise(mean = mean(x),
lower = quantile(x, 0.025),
upper = quantile(x, 0.975),
.groups = "drop") %>%
mutate(time = target_data_KONZ_sub$time)
plot_TempDLM_KONZ <- ggplot(data = plot_chain_TempDLM_KONZ, aes(x = time, y = mean)) +
geom_line() +
geom_ribbon(aes(ymin = lower, ymax = upper),
alpha = 0.2,
color = plot_colors[2],
fill = plot_colors[2]) +
geom_point(data = target_data_KONZ_sub, aes(x = time, y = nee), color = "black") +
labs(title = paste("Site: KONZ | Time Period: ", substr(target_data_KONZ_sub$time[1], 1, 10), "to", substr(target_data_KONZ_sub$time[subset_length_timesteps], 1, 10)),
x = "Time (UTC)",
y = expression(paste("NEE (",mu, "mol" ~ CO[2] ~ m^-2 ~ s^-1,")"))) +
theme_bw()
# OSBS
plot_chain_TempDLM_OSBS <- chain_TempDLM_OSBS %>%
group_by(day) %>%
summarise(mean = mean(x),
lower = quantile(x, 0.025),
upper = quantile(x, 0.975),
.groups = "drop") %>%
mutate(time = target_data_OSBS_sub$time)
plot_TempDLM_OSBS <- ggplot(data = plot_chain_TempDLM_OSBS, aes(x = time, y = mean)) +
geom_line() +
geom_ribbon(aes(ymin = lower, ymax = upper),
alpha = 0.2,
color = plot_colors[3],
fill = plot_colors[3]) +
geom_point(data = target_data_OSBS_sub, aes(x = time, y = nee), color = "black") +
labs(title = paste("Site: OSBS | Time Period: ", substr(target_data_OSBS_sub$time[1], 1, 10), "to", substr(target_data_OSBS_sub$time[subset_length_timesteps], 1, 10)),
x = "Time (UTC)",
y = expression(paste("NEE (",mu, "mol" ~ CO[2] ~ m^-2 ~ s^-1,")"))) +
theme_bw()
# SRER
plot_chain_TempDLM_SRER <- chain_TempDLM_SRER %>%
group_by(day) %>%
summarise(mean = mean(x),
lower = quantile(x, 0.025),
upper = quantile(x, 0.975),
.groups = "drop") %>%
mutate(time = target_data_SRER_sub$time)
plot_TempDLM_SRER <- ggplot(data = plot_chain_TempDLM_SRER, aes(x = time, y = mean)) +
geom_line() +
geom_ribbon(aes(ymin = lower, ymax = upper),
alpha = 0.2,
color = plot_colors[4],
fill = plot_colors[4]) +
geom_point(data = target_data_SRER_sub, aes(x = time, y = nee), color = "black") +
labs(title = paste("Site: SRER | Time Period: ", substr(target_data_SRER_sub$time[1], 1, 10), "to", substr(target_data_SRER_sub$time[subset_length_timesteps], 1, 10)),
x = "Time (UTC)",
y = expression(paste("NEE (",mu, "mol" ~ CO[2] ~ m^-2 ~ s^-1,")"))) +
theme_bw()
plot_TempDLM_BART
plot_TempDLM_KONZ
plot_TempDLM_OSBS
plot_TempDLM_SRER
Compare temperature DLM to random walk fit. As can be seen, the two models only slightly differ when there is a relatively large gap in the data.
# BART
comparison_plot_BART <- ggplot() +
geom_line(data = plot_chain_RW_BART, aes(x = time, y = mean, color = "1")) +
geom_line(data = plot_chain_TempDLM_BART, aes(x = time, y = mean, color = "2")) +
geom_ribbon(data = plot_chain_RW_BART, aes(x = time, ymin = lower, ymax = upper, color = "1", fill = "1"), alpha = 0.2) +
geom_ribbon(data = plot_chain_TempDLM_BART, aes(x = time, ymin = lower, ymax = upper, color = "2", fill = "2"), alpha = 0.2) +
geom_point(data = target_data_BART_sub, aes(x = time, y = nee), color = "black") +
scale_fill_manual(name = "Model", values = c("gray", plot_colors[1]), labels = c("Random Walk", "Temp. DLM")) +
scale_color_manual(name = "Model", values = c("gray", plot_colors[1]), labels = c("Random Walk", "Temp. DLM")) +
labs(title = paste("Site: BART | Time Period: ", substr(target_data_BART_sub$time[1], 1, 10), "to", substr(target_data_BART_sub$time[subset_length_timesteps], 1, 10)),
x = "Time (UTC)",
y = expression(paste("NEE (",mu, "mol" ~ CO[2] ~ m^-2 ~ s^-1,")"))) +
theme_bw()
# KONZ
comparison_plot_KONZ <- ggplot() +
geom_line(data = plot_chain_RW_KONZ, aes(x = time, y = mean, color = "1")) +
geom_line(data = plot_chain_TempDLM_KONZ, aes(x = time, y = mean, color = "2")) +
geom_ribbon(data = plot_chain_RW_KONZ, aes(x = time, ymin = lower, ymax = upper, color = "1", fill = "1"), alpha = 0.2) +
geom_ribbon(data = plot_chain_TempDLM_KONZ, aes(x = time, ymin = lower, ymax = upper, color = "2", fill = "2"), alpha = 0.2) +
geom_point(data = target_data_KONZ_sub, aes(x = time, y = nee), color = "black") +
scale_fill_manual(name = "Model", values = c("gray", plot_colors[2]), labels = c("Random Walk", "Temp. DLM")) +
scale_color_manual(name = "Model", values = c("gray", plot_colors[2]), labels = c("Random Walk", "Temp. DLM")) +
labs(title = paste("Site: KONZ | Time Period: ", substr(target_data_KONZ_sub$time[1], 1, 10), "to", substr(target_data_KONZ_sub$time[subset_length_timesteps], 1, 10)),
x = "Time (UTC)",
y = expression(paste("NEE (",mu, "mol" ~ CO[2] ~ m^-2 ~ s^-1,")"))) +
theme_bw()
# OSBS
comparison_plot_OSBS <- ggplot() +
geom_line(data = plot_chain_RW_OSBS, aes(x = time, y = mean, color = "1")) +
geom_line(data = plot_chain_TempDLM_OSBS, aes(x = time, y = mean, color = "2")) +
geom_ribbon(data = plot_chain_RW_OSBS, aes(x = time, ymin = lower, ymax = upper, color = "1", fill = "1"), alpha = 0.2) +
geom_ribbon(data = plot_chain_TempDLM_OSBS, aes(x = time, ymin = lower, ymax = upper, color = "2", fill = "2"), alpha = 0.2) +
geom_point(data = target_data_OSBS_sub, aes(x = time, y = nee), color = "black") +
scale_fill_manual(name = "Model", values = c("gray", plot_colors[3]), labels = c("Random Walk", "Temp. DLM")) +
scale_color_manual(name = "Model", values = c("gray", plot_colors[3]), labels = c("Random Walk", "Temp. DLM")) +
labs(title = paste("Site: OSBS | Time Period: ", substr(target_data_OSBS_sub$time[1], 1, 10), "to", substr(target_data_OSBS_sub$time[subset_length_timesteps], 1, 10)),
x = "Time (UTC)",
y = expression(paste("NEE (",mu, "mol" ~ CO[2] ~ m^-2 ~ s^-1,")"))) +
theme_bw()
# SRER
comparison_plot_SRER <- ggplot() +
geom_line(data = plot_chain_RW_SRER, aes(x = time, y = mean, color = "1")) +
geom_line(data = plot_chain_TempDLM_SRER, aes(x = time, y = mean, color = "2")) +
geom_ribbon(data = plot_chain_RW_SRER, aes(x = time, ymin = lower, ymax = upper, color = "1", fill = "1"), alpha = 0.2) +
geom_ribbon(data = plot_chain_TempDLM_SRER, aes(x = time, ymin = lower, ymax = upper, color = "2", fill = "2"), alpha = 0.2) +
geom_point(data = target_data_SRER_sub, aes(x = time, y = nee), color = "black") +
scale_fill_manual(name = "Model", values = c("gray", plot_colors[4]), labels = c("Random Walk", "Temp. DLM")) +
scale_color_manual(name = "Model", values = c("gray", plot_colors[4]), labels = c("Random Walk", "Temp. DLM")) +
labs(title = paste("Site: SRER | Time Period: ", substr(target_data_SRER_sub$time[1], 1, 10), "to", substr(target_data_SRER_sub$time[subset_length_timesteps], 1, 10)),
x = "Time (UTC)",
y = expression(paste("NEE (",mu, "mol" ~ CO[2] ~ m^-2 ~ s^-1,")"))) +
theme_bw()
comparison_plot_BART
comparison_plot_KONZ
comparison_plot_OSBS
comparison_plot_SRER